# -*- coding: utf-8 -*-
"""
Created on Fri May 31 20:22:52 2013

@author: Maxim
"""

import airfoil
from numpy import array,zeros,ones
from ga2 import *
import ga
from scipy.interpolate import Rbf,interp1d
from scipy.optimize import fminbound
import matplotlib.pyplot as plt


class Design:
    def __init__(self):
        self.af = airfoil.Airfoil()
        self.Mcrs = 0.179383051018
        self.Mldg = 0.0940384447251
        self.Recrs = 4415299.0
        self.Reldg = 2700149.0
        self.alphaSeq = [-5.0,20.0,1.0]
        self.clreq = [0.2,0.4,0.6]
        self.clmaxReq = 1.3
        self.tcMin = 0.14
        self.tcMax = 0.16
        self.alphaClmax = 15.0
        self.boattailMin = 10.0
    def set_cst(self,x):
        Au = array([x[0],x[1],x[2],x[3]])
        Al = array([-x[0],x[4],x[5],x[6]])
        self.af.create_CST(Au,Al)
        self.boattail = ny.degrees(ny.arctan(Au[-1])-ny.arctan(Al[-1]))
    def objective(self,x,params=None):
        g = zeros(5)
        try:
            self.set_cst(x)
            polar1 = self.af.get_X_polar(self.Mcrs,self.Recrs,self.alphaSeq,50)
            polar1.create_splines()
            cd = zeros(len(self.clreq))
            for i,cl in enumerate(self.clreq):
                cd[i] = polar1.get_cd_at_cl(cl)
            polar2 = self.af.get_X_polar(self.Mldg,self.Reldg,self.alphaSeq,50)
            polar2.create_splines()
            polar2.calc_clmax()
            g[0] = self.clmaxReq - polar2.clmax
            g[1] = self.tcMin - self.af.thickness
            g[2] = self.af.thickness - self.tcMax
            g[3] = self.alphaClmax - polar2.alphaClmax
            g[4] = self.boattailMin - self.boattail
            f = cd.mean()
            return f, g,0.0
        except:
            return 1e5,-100.0*ones(len(g)),0.0
    def finalize(self,x):
        self.set_cst(x)
        self.af.write_txt('airfoil_result_09.txt')
        self.af.plot()
        pol = self.af.get_X_polar(self.Mldg,self.Reldg,self.alphaSeq,50)
        pol.display()

def stage_1():
    lb = array([0.08,0.001,0.001,0.001,-0.4   ,-0.4   ,-0.4])
    ub = array([0.3 ,0.4  ,0.4  ,0.4  , 0.1, 0.1, 0.1])
    design = Design()
    cost = design.objective
    
    opt = gaOptions(lb,ub)
    opt.CorrNan = True
    opt.PopSize = 100
    opt.MaxIterations(150,75)
    opt.sigmaEnd = .1
    opt.sigma = 0.9
    opt.MutRate = 0.2
    opt.EliteRatio = 0.2
    opt.HistFile = 'airfoil_history_09.txt'
    fBest,xBest,fHist,xHist,Iter = gaMain(cost, opt)
    design.finalize(xBest)

class TmpPolar:
    def __init__(self):
        self.alpha = list()
        self.cl = list()
        self.cd = list()
        self.cm = list()
    def add(self,alpha,cl,cd,cm):
        self.alpha.append(alpha)
        self.cl.append(cl)
        self.cd.append(cd)
        self.cm.append(cm)
    def _create_splines(self):
        self.lift = interp1d(self.alpha,self.cl,'cubic')
        self.drag = interp1d(self.alpha,self.cd,'cubic')
        self.moment = interp1d(self.alpha,self.cm,'cubic')
    def get_clmax(self):
        self._create_splines()
        def obj(alpha):
            return -self.lift(alpha)
        self.alphaClmax = fminbound(obj,self.alpha[1], self.alpha[-2])
        self.clmax = self.lift(self.alphaClmax)
        return self.clmax
    def plot(self):
        plt.figure(1)
        plt.grid(True)
        plt.plot(self.alpha,self.cl,'bo-')
        plt.ylabel('lift coefficient')
        plt.xlabel('alpha,deg')
        plt.figure(2)
        plt.grid(True)
        plt.plot(self.cd,self.cl,'ro-')
        plt.xlabel('drag coefficient')
        plt.ylabel('lift coefficient')
#        plt.figure(3)
#        plt.grid(True)
#        plt.plot(self.alpha,self.cm,'go-')
#        plt.xlabel('alpha,deg')
#        plt.ylabel('moment coefficient')
        plt.show()
        
def read_polar(path):
    fid = open(path,'rt')
    for i in range(4):
        fid.readline()
    pol = TmpPolar()
    for line in fid:
        if line.strip()!='':
            seg = line.split()
            seg = [float(s) for s in seg]
            pol.add(seg[0],seg[1],seg[2],seg[3])
    return pol

def normalize(x,bnds):
    xnorm = ny.zeros(len(x))
    for i,xx in enumerate(x):
        lb = bnds[i][0]
        ub = bnds[i][1]
        xnorm[i] = (2.0*xx-(lb+ub))/(ub-lb)
    return xnorm

def denormalize(xnorm,bnds):
    x = ny.zeros(len(xnorm))
    for i,xx in enumerate(xnorm):
        lb = bnds[i][0]
        ub = bnds[i][1]
        x[i] = 0.5*(xnorm[i]*(ub-lb)+(lb+ub))
    return x

class Design2:
    def __init__(self,xDoe=None,clmaxDoe=None):
        if xDoe!=None and clmaxDoe!=None:
            self.clmaxModel = Rbf(xDoe[:,0],xDoe[:,1],xDoe[:,2],xDoe[:,3],xDoe[:,4],clmaxDoe)
        self.af = airfoil.Airfoil()
        self.Mcrs = 0.179383051018
        self.Recrs = 4415299.0
        self.alphaSeq = [-10.0,20.0,1.0]
        self.clreq = [0.2,0.4,0.6]
        self.cdMin = [0.004270,0.004327,0.004955]
        self.tcMin = 0.12
        self.tcMax = 0.16
        self.Au0 = 0.185091
        self.Al3 = 0.059345
    def set_cst(self,x):
        Au = array([self.Au0,x[0],x[1],x[2]])
        Al = array([-self.Au0,x[3],x[4],self.Al3])
        self.af.create_CST(Au,Al)
    def _tmp_run(self,x):
        self.set_cst(x)
        polar1 = self.af.get_X_polar(self.Mcrs,self.Recrs,self.alphaSeq,50)
        polar1.create_splines()
        cd = zeros(len(self.clreq))
        for i,cl in enumerate(self.clreq):
            cd[i] = polar1.get_cd_at_cl(cl)
        print cd
    def objective(self,x):
        try:
            self.set_cst(x)
            polar1 = self.af.get_X_polar(self.Mcrs,self.Recrs,self.alphaSeq,50)
            polar1.create_splines()
            cd = zeros(len(self.clreq))
            for i,cl in enumerate(self.clreq):
                cd[i] = polar1.get_cd_at_cl(cl)
            g = zeros(3)
            polar2 = self.af.get_X_polar(self.Mldg,self.Reldg,self.alphaSeq,50)
            polar2.create_splines()
            polar2.calc_clmax()
            g[0] = self.clmaxReq - polar2.clmax
            g[1] = self.tcMin - self.af.thickness
            g[2] = self.af.thickness - self.tcMax
            return cd.sum(),g*100.0,0.0
        except:
            return 100.0,-100.0*ones(3),0.0
    def objective2(self,x):
        try:
            self.set_cst(x)
            polar = self.af.get_X_polar(self.Mcrs,self.Recrs,self.alphaSeq,50)
            polar.create_splines()
            cd = zeros(len(self.clreq))
            for i, cl in enumerate(self.clreq):
                cd[i] = polar.get_cd_at_cl(cl)
            g = zeros(5)
            g[:3] = (cd - self.cdMin)*1000.0
            g[3] = self.tcMin - self.af.thickness
            g[4] = self.af.thickness - self.tcMax
            return -self.clmaxModel(x[0],x[1],x[2],x[3],x[4]),g,0.0
        except:
            return 1000.0,-1000*ones(len()),0.0

    def finalize(self,x):
        self.set_cst(x)
        self.af.write_txt('airfoil_result_05.txt')
        self.af.plot()

def postprocess_2():
    x = array([0.259564,0.289934,0.195339,-0.144149,-0.232191])
    des = Design2()
    des._tmp_run(x)

def stage_2():
    clmaxPath = r'D:\4. Conferences\AIAA CFD 2013\af_results\DoE_CCD'
    DoePath = r'D:\3. Projects\afOpt\source\DoE_desings_20130603.txt'
    x0 = array([0.257832,0.274773,0.187216,-0.155010,-0.239510])
    dx = 0.02
    bnds = ny.transpose(ny.vstack([x0-dx,x0+dx]))
    nDoE = 44
    xDoE = ny.zeros([nDoE,5])
    clmaxDoE = ny.zeros(nDoE)
    for i in range(nDoE):
        pol = read_polar(clmaxPath+'\\airfoil_doe_%d.txt'%(i+1))
        clmaxDoE[i] = pol.get_clmax()
    fid = open(DoePath,'rt')
    fid.readline()
    for i in range(nDoE-1):
        line = fid.readline()
        seg = line.split()
        x = [float(s) for s in seg]
        xDoE[i] = denormalize(x,bnds)
    xDoE[-1] = array([0.261066,0.291011,0.197810,-0.143179,-0.232682])
    
    def equal(arr1,arr2):
        for a1,a2 in zip(arr1,arr2):
            if a1!=a2:
                return False
        else:
            return True
    
    for i,xx in enumerate(xDoE):
        for j,x2 in enumerate(xDoE):
            if i!=j and equal(xx,x2):
                print '%d is equal to %d'%(i,j)
                print clmaxDoE[i], clmaxDoE[j]
    
    design = Design2(xDoE,clmaxDoE)
    cost = design.objective2
    
    opt = gaOptions(x0-dx,x0+dx)
    opt.CorrNan = True
    opt.PopSize = 75
    opt.MaxIterations(100,100)
    opt.sigmaEnd = .1
    opt.sigma = 1.0
    opt.HistFile = 'airfoil_history_05.txt'
    fBest,xBest,fHist,xHist,Iter = gaMain(cost, opt)

    design.finalize(xBest)

def coord_cst(x):
    au = x[:4]
    al = array([-x[0],x[4],x[5],x[6]])
    af = airfoil.Airfoil()
    af.create_CST(au,al)
    return af.coord

def plot_results():
    af1 = airfoil.Airfoil()
    af2 = airfoil.Airfoil()
    af3 = airfoil.Airfoil()
    af3.read_txt('GA37A315mod_reworked.txt')
    Mcrs = 0.179383051018
    Recrs = 4415299.0
    x1 = array([0.185091,0.257832,0.274773,0.187216,-0.155010,-0.239510,0.059345])
    x2 = array([0.185091,0.259564,0.289934,0.195339,-0.144149,-0.232191,0.059345])
    af1.create_CST(x1[:4],[-x1[0],x1[4],x1[5],x1[6]])
    af2.create_CST(x2[:4],[-x2[0],x2[4],x2[5],x2[6]])
    pol1 = af1.get_X_polar(Mcrs,Recrs,[-10,20,0.5],50)
    pol2 = af2.get_X_polar(Mcrs,Recrs,[-10,20,0.5],50)
    pol3 = af3.get_X_polar(Mcrs,Recrs,[-10,20,0.5],50)
    plt.figure(1)
    plt.grid(True)
    plt.hold(True)
    plt.axis('equal')
    plt.plot(af1.coord[:,0],af1.coord[:,1],'b-')
    plt.plot(af2.coord[:,0],af2.coord[:,1],'r-')
    plt.legend(['1st stage optimum','2nd stage optimum'])

    plt.figure(2)
    plt.grid(True)
    plt.hold(True)
    plt.plot(pol1.cd,pol1.cl,'b-')
    plt.plot(pol2.cd,pol2.cl,'r-')
    #plt.plot(pol3.cd,pol3.cl,'g-')
    plt.legend(['1st stage optimum','2nd stage optimum'],'lower right')
    plt.xlabel('drag coefficient')
    plt.ylabel('lift coefficient')
    plt.xlim([0.0,0.03])
    plt.show()

def stage_11():
    lb = array([0.05,0.001,0.001,0.001,-0.4   ,-0.4   ,-0.4])
    ub = array([0.3 ,0.4  ,0.4  ,0.4  , 0.1, 0.1, 0.1])
    design = Design()
    costFcn = design.objective
    GA = ga.ga(costFcn,lb,ub)
    GA.setCPU(1)
    GA.displayFlag = 1
    GA.solve()
    print 'f=%.4f x1=%.4f x2=%.4f, Neval=%d' %(GA.fOpt,GA.xOpt[0],GA.xOpt[1],GA.funcCount)

if __name__=="__main__":
    stage_1()
